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Abstract 

Generalized Estimation Equations (GEE) are a well-known method for the anal¬ 
ysis of non-Gaussian longitudinal data. This method has computational simplic¬ 
ity and marginal parameter interpretation. However, in the presence of missing 
data, it is only valid under the strong assumption of missing completely at random 
(MCAR). Some corrections can be done when the missing data mechanism is missing 
at random (MAR): inverse probability weighting (WGEE) and multiple imputation 
(MIGEE). In order to obtain consistent estimates, it is necessary the correct speci¬ 
fication of the weight model for WGEE or the imputation model for the MIGEE. A 
recent method combining ideas of these two approaches has doubly robust property. 

For consistency, it requires only the weight or the imputation model to be correct. In 
this work it is assumed a proportional odds model and it is proposed a doubly robust 
estimator for the analysis of ordinal longitudinal data with intermittently missing 
response and covariate under the MAR mechanism. Simulation results revealed 
better performance of the proposed method compared to WGEE and MIGEE. The 
method is applied to a data set related to Analgesia Pain in Childbirth study. 

Keywords: Missing at random; Multiple imputation; Proportional Odds Model; Weighted 
GEE. 


1 Introduction 


The Generalized Estimating Equation (GEE) method (Liang & Zeger 1986) is one 
of the most popular approaches for the analysis of non-Gaussian correlated data. Its 
main advantage resides in the fact that one is only required to specify correctly the mean 
structure of the response for the parameter estimator to be consistent and asymptotically 
normal. In its basic formulation the association parameters among repeated measures 
were taken as nuisance parameters. GEE method has computational simplicity (there is 
no need of dealing with complex, and in some cases, intractable likelihoods) and it further 
allows populational-averaged interpretation of the parameter of interest. 
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It is very common for sets of longitudinal data to be incomplete, in the sense that not 
all planned observations are actually observed. This problem is pervasive in longitudinal 
data because nonresponse can occur any time from the beginning of the study. Two 
patterns of missing data can be observed for the response: (1) dropout, when a subject 
leaves the study prematurely for reasons beyond the control of the investigator, leading to 
a monotone pattern of nonresponse, or (2) intermittent nonresponse, in which a subject 
returns to the study after some occasions of nonresponse. Covariates may also be missing, 
leading to limitations in data analysis. In the presence of missing data, three issues are 
of main concern: (1) potential serious bias due to systematic differences between the 
observed data and the missing data, (2) complications in data handling and statistical 
inferences, and (3) loss of efficiency. Therefore, in order to make valid inferences it is 
fundamental to know the missing data mechanism generating the nonresponse and how 
to handle it. 


Little & Rubin (1987) provided a formal framework for dealing with missing data by 


defining the commonly adopted taxonomy of missing data mechanisms. A nonresponse 
process is said to be missing completely at random (MCAR) if missingness is independent 
of both unobserved and observed data, and missing at random (MAR) if, conditionally 
on the observed data, the missingness is independent of the unobserved data. When the 
nonresponse process depends on unobserved quantities it is said to be missing not at 
random (MNAR). 

When data are incomplete, GEE suffers from its frequentist nature and is, in its basic 
form, valid only under MCAR (Liang & Zeger, 1986). The first effort to make GEE appli¬ 
cable to the more realistic MAR scenario was Multiple Imputation (MIGEE), proposed by 


Little & Rubin (1987), in which the missing portions of data are multiply imputed taking 
into account the uncertainty associated with the predicted values. The completed data 
sets are analyzed by standard methods for complete data, and estimates are combined into 
a final analysis. Multiple imputation is detailed in the books by Schafer (|1997), Little & 


Rubin (2002) and Carpenter & Kenward (2013). Later, Robins et al. (1995) proposed the 


Weighted Generalized Estimating Equations (WGEE), which consists in weighting each 
observation by the inverse of the probability of the data being observed. This method 
produces consistent estimates provided the weight model is correctly specified. 

Doubly robust estimators (DRGEE) arise as a third generalization of ordinary GEE 
to deal with data subject to MAR mechanism. DR methods have received increasingly 


attention in the literature in the last decade (see Carpenter et al. (2006), Bang & Robins 


(2005), Tsiatis (2006), Seaman & Copas (2009), Chen & Zhou (2011)). The main idea is 
to supplement the WGEE with a predictive model for the missing quantities conditional 
on the observed ones. Doubly Robust method requires only the dropout or the conditional 
model to be correctly specified in order to provide consistent estimates. In the analysis 
of longitudinal binary data, DR methods have been applied by Seaman & Copas (2009), 


Birhanu et al. (2011), for missing responses and by Chen & Zhou (2011) for intermittently 


missing response and a single missing covariate. 

Literature of GEE for missing data is comparatively scarce for longitudinal ordinal 
response. In Toledano & Gatsonis (1999), the authors used a weighted GEE method to 
accommodate arbitrary patterns of a MCAR missing response and missingness in a key 
covariate subject to a MAR mechanism. A recent paper from Donneau et al. (2014a) com¬ 
pared through a simulation study two multiple imputation methods (multivariate normal 
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imputation and ordinal imputation regression) for longitudinal ordinal data subject to 
dropout. In another paper the same authors compared joint modeling and fully condi¬ 
tional specification approaches for non-monotone missingness (Donneau et al., 2014b). 


The above mentioned papers used single robust versions of GEE and they have treated 
only a missing MAR baseline covariate or missing MAR response. Thus the use of a DR 
GEE method for ordinal data with simultaneously intermittently missing response and 
missing covariate has been in need of further development. 

This work was motivated by the Analgesia in Childbirth study which was conducted 
in Minas Gerais state, Brazil. The main objective of that study was to compare two 
techniques of analgesia for labor pain in 49 patients. The response, pain intensity, was 
subjectively assessed by each patient, and various clinical covariates were taken until 
delivery. Response and a particular covariate (consumption of oxytocin) were missing for 
some patients and the MAR mechanism seems to be a reasonable assumption for this 
data. 

In the current paper, it is proposed a doubly robust approach for the analysis of 
longitudinal ordinal data with intermittently missing response and covariate that is MAR. 
The proposed methodology, to our knowledge new to the GEE literature of missing ordinal 
data, can be used for handling arbitrary patterns of missing data in the ordinal response 
and missingness in a key covariate, as those frequently arising in medical studies. 

The paper is organized as follows. In Section [2] are defined the notation for GEE with 
fully observed data and missing data mechanisms. Section [3] outlines WGEE and MIGEE 
approaches. The proposed methodology is established in Section [4j A simulation study is 
presented in Section [5j in which the finite-sample biases and standard errors are compared 
for the standard GEE, MIGEE, WGEE and doubly robust versions. Data arising from 
the Analgesia en Childbirth study are analyzed in Section [6j Paper ends with a discussion 
and future directions in Section 0 


2 GEE for Complete Data and Missing Data As¬ 
sumptions 


In this section it is introduced generalized estimating equations for the analysis of fully 
observed ordinal data. Section M establishes the model and notation for longitudinal 
ordinal data. Section A2 presents a series of assumptions related to mechanism causing 
data to be missing and necessary to be considered in order to build valid estimators. 


2.1 GEE for Longitudinal Ordinal Response 

Let On G {1, 2,..., J} be the ordinal response for subject i (i = 1,..., n) at time 
t (t = 1,..., Tj, Ti <T). As the response has J levels it can be defined Y itj = I {On = j) 
for j = 1,..., J, where I {A) denotes the indicator function. Y it] is converted into the 
equivalent (J — l)-variate vector Y it = (Y in ,..., Y it (j_ p) r and let Y = (Y^f,..., Y^J T 
the stacked response vector. When J = 2 the response is binary and Y it is a scalar. Let 
Xi = (A"^,..., Xj T .) T denotes the T t x 1 covariate vector that may be missing for the i-th 
subject, and Z r = (Zj l7 ..., Zj Tj ) T the Ti x q matrix of explanatory variables that are 
always observed. 
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The marginal distribution of Y it is assumed to be multinomial (with sample size 
J2j =1 Y itj = 1), that is 

j 

f(Y it \X it ,Z lt ,f3) = (1) 

3 =1 

where p itj = puj((3) = E(Y itj \X il Z h (3) = Pr(O it = j\Xi,Zi,[3), is the probability of 
response j at time t and /3 is a p x 1 vector of parameters. Two common choices for 
modeling fi it j are the cumulative logit and probit models. In this work it is assumed a 
cumulative logit link, that is, 


lo&t[Pr(Oit<j\X it ,Z it )]=p Qj + X it p x + ZT t /3 g , j = (2) 

Formulation in ([2]) implies a proportional odds model (McCullagh, 1980). In such 
model the interpretation of (3 is the same regardless of the number of categories (i.e., it is 
invariant to combination of categories). A desired feature is that the exponential of the 
parameters is interpreted as an odds ratio (Agresti, 2013). 

Main interest is to make inferences related to the regression parameters 
(3 = (/?oi, • • •, f3o,j-i, P x , /3T) t associated with the (J — 1) x 1 marginal probability vectors 

E{Yn\X i, Z j) /Xj t (/3) (/^ith ■ ■ • j Mit(j-i)) ■ 

fj, it is grouped to form a vector E(Yi\Xi, Z{) — pb t — (/x^,..., fjJ Ti ) T with the same 
dimension of Yi. 

In order to estimate /3 generalized estimation equations are used (Liang & Zeger (1986); 


Lipsitz et al. (1994)), which takes the form 

n n 

U(I3) = J2 um = d,v-\y, - m.) = 0, 


(3) 


i =1 


2—1 


where D % = 0- and Vi = Vi(/3, a) is a T)(J — 1) x T)(J — 1) “working covariance” 

matrix usually decomposed into the form Vi(/3,ot ) = F 1 / 2 ((3)C{(a.)F 1 / 2 ((3) , where F, is 
a matrix containing the marginal variances, F it , given by 


Fit diag \_Pit 1 (1 pitl );•••; Pit,J— 1 (1 Pit, J— l)] , 

and Ci is equal to the marginal correlation matrix. The (J — 1) x (J — 1) diagonal blocks 
of Ci(a ) are F~ t 1/2 V it F~ t 1/2 , with V it = diag(/x it ) - and the ( J - 1) x ( J - 1) 

off-diagonal blocks of (7* (a) are am/, which represents the correlation between Y lt and 
Ya', t^t' (Lipsitz et al, |1994 ). 

Under mild regularity conditions and correct specification of the marginal mean model 
in (J2]) , Liang and Zeger (1986) proved that the estimator /3, obtained by solving (jij), is 
consistent and y/n((3 — /3) converges in distribution to a p -variate normal distribution with 
mean 0 and covariance matrix 


(Lipsitz et al., |199 4 


V p = lim uS 0 1 SiS f 


(4) 


where E r 


= V n DV l D T 

L^i= 1 ' v i ’ 


and Ex = Y^i=\ ^iVi 1 Cov(Y i )‘\Z i 1 Dj . In practice, the 


“sandwich” covariance matrix in Q is calculated by ignoring the limit and replacing 
((3, a) and Cov(lU) by (/ 3, a ) and (Y i — /2 i )(Y' \ — /2j) r , respectively (Touloumis et al 


2013). 


4 


















2.2 Missing Data Framework 

For each occasion t it can be defined Ru = 0 if Ou and X it are missing, R lt = 1 if O lt 
is missing and X it is observed, Ru = 2 if Ou is observed and X lt is missing, and Ru = 3 
if O it and X lt are both observed. Let R, = (Ru, ..., R^t,• i ) T , and R lt = (Ru ,..., R^t- 1 )- 

By specifying conditional models of the form Pr(R it = r it (),. XZ,) it can be 
obtained Pr(R t = r^O^X,, Z { ) through ]l2=2 jPr (- R « = ru\Ru,Oi, X { , Z i )Pr(R il = 
rn\Oi,Xi,Zi). Let \ ltk = Pr(R it = k\R it ,Oi, X i} Z { ), for k = 0,1,2,3. This general 
formulation encompasses MCAR, MAR and MNAR mechanisms. In particular, the MAR 
mechanism requires 


Pr(R t 


r,\O u X h Zi) = Pr(Ri = n\0°, X°, Z lt ), 


(5) 


where 0° and X° denotes the observed components of O l and X,, respectively. Because 
R, is modeled through a product of conditional models, it is natural to make the following 


further assumption (Chen & Zhou, 2011) 


Pr(R it = r it \Ru, O b , X h Z { ) = Pr(R it = r it \Ru, 0° it , X° it , Z % 


0 , 


( 6 ) 


for each time t, where 0° it and X° it are the histories of observed responses and covariates 
up to time t — 1. 

Let itu = Pr(Ru = 3|Oj, X t , Zi) be the marginal probability of observing both Oi and 
Xi at time t, given the entire vectors of responses and covariates. Then, 7i it is expressed 
by 

^it ^ 3, Ri i— i , Ru rn\Oi, Xi, Z ^). 

This marginal probability can be expressed in terms of the conditional probabilities A^’s. 
Throughout the paper it is required the so-called positivity assumption, that is, 7 T it must 
be bounded away from zero. This condition is needed in order to guarantee the existence 
of -y/n-consistent estimators of [3 


(Robins et al ., 1995 


3 Available Approaches for Missing Data 


Multiple imputation and weighted generalized estimation equations are two commonly 
used methods available for missing data under MAR mechanism. These methods are 


presented in Sections 3T and T2, respectively. They serve as the basis for the construction 
of the doubly robust estimator, presented in Section |4j 


3.1 Multiple Imputation Generalized Estimating Equations 

A imputation model commonly used to handle intermittently missing response and 
covariate, is imputation using chained equations (van Buuren et al. ( 1999[ ), van Buuren 
(2007)), which is more commonly referred to as full conditional specification (FCS). This 


approach specifies conditional distributions for each incomplete variable, conditional on 
all others variables in the imputation model. Starting from an initial imputation, FCS 
draws imputations by iterating over the conditional densities. 
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Denote by (3 m and U m , respectively, the estimate of (3 and its covariance matrix from 


the GEE analysis of the m-th completed data set, (m = 1, Following Rubin 


(1987), the combined point estimate for the parameter of interest /3 from the MI is simply 


the average of the M complete-data point estimates 


1 


M 


PMI ^ ^ fimi 

m =1 


and an estimate of the covariance matrix of (3 MI is given by 

' M + 1' 


where 


1 M 

m =1 


U Ml = w 


and B = 


M 


B, 


M 


M - 1 


^2 (Pm ~ &Ml)((3 m ~ Pmi)' ■ 


m =1 


3.2 Weighted Generalized Estimating Equations 


Robins et al. ([1995 ) proposed a class of weighted estimating equations to allow for 


MAR mechanism. In binary longitudinal data, Chen & Zhou (2011) extended the method 


to accommodate arbitrary patterns of missing response and missing covariate. Their 
method was adapted here for longitudinal ordinal responses. 

Define a weight matrix A* = t = 1, • • •, T h t! = 1 ,..., T;, where 

<W = {I{Ru = l, Rw = 3) + I(Ri t = 3 ,Rif = 3 )}/ 7 T itt ' for t 7 ^ t', 8 itt = I(R it = 3)/n it , 
and 7 Tut' = Pr(R it = 1 ,R it > = 3| O u X h Zi) + Pr(R it = 3 ,R it > = 3\O t , X t , Zj). Let 
Mi = F i 1//2 (C“ 1 • A ijFl 1,72 where A • B = [a it ■ b lt ] denotes the Hadamard product of 
matrix A = [an] and B = [bn]. 

The weighted generalized estimating equations (WGEE) for /3 are given by 


c/(/3,^) = ^c/ i (A^) = 0, 


(7) 


i= 1 


where U i ([3,'ip) = D M,{Y t - /zj. A consistent estimate for /3 can be obtained by 
solving (Jt[) , under the correct specification of the missing data model. 

To model \ it k it is adopted a politomic logistic regression, with Aito as the reference 
category, that is 

/Xi±A - - - (8) 


log 1 W ) = u itk^ki k = 1,2,3, 

V Mo J 

where the covariates u itk are some function of {Ri t ,0° t ,X° t , Z it }. 


4 Doubly Robust GEE for Longitudinal Ordinal Data 


Some authors (e.g., Scharfstein et al. ( 1999[ ) , |Tsiatis| ( [20 06 )) noted that adding a term 
of expectation zero, say </>(•), to the inverse probability weighted estimators would still 
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result in consistent estimates under a MAR mechanism. The solutions of these augmented 
estimating equations give rise to the so-called doubly robust estimators. 


Chen & Zhou (2011) showed that the optimal (j) op t for missing response and covariate is 

given by (p opt = E(YY L ,x r i n \Y°,x?,z i ,R i ) {DiN^Yi — /zj}, with AC = F ■ 1//2 [C i 1 • (11 T — A*)} F ■ 1//2 , 
where 1 is a vector of l’s of length Tj(J — 1), and Y™ and A” 1 denote the missing com¬ 
ponents of Yi and X i} respectively. 

An improved estimate for /3 can then be obtained by solving the estimating equations 

n n 

Si{0) = = - /zj + E(Y™ ) x™\Y°,x°,z,,R l ) - Hi)}] = 0. 




i =1 


(9) 

The estimator for / 3 in (J9]) is doubly-robust in the sense that it is consistent if at least one 
of the missing data model or the covariate model is correctly specified. 

Applications of doubly robust estimators in longitudinal settings include Bang & 


Robins (2005), Seaman & Copas (2009), Chen & Zhou (2011) and Birhanu et al. ( 2011[ ). 
Those developments focus mainly on binary response and have not been, to our knowledge, 
investigated with ordinal longitudinal data. In this work it is considered a longitudinal 
response measured on a ordinal scale. 

The referred expectation in the second part of (|9]) is over the conditional distribution 
of ( Y - n , X™\Y°, X°, Zi, Rj ), which can be written as 

P(YT = V?:X? = x?\Y°,X°,Z, n R fl (3* n ) = P(Y? = y?,X? = x?\ Y^X^Z*?^) 

= P(Y? = yr\Y° i ,X i = x i ,Z i ;F) 
xP{X? = x?\Y° i ,X° i ,Z i] 7 ). 

The multivariate distribution P(Y™ = y™\Y°, Xi = x^ Zi,(3*) is expressed through 
a product of univariate ordinal models. 


4.1 Estimation for the Nuisance Parameters 

The method of maximum likelihood is employed to estimate t/u The log likelihood of 
the politomic logistic model for i/? has the form 

n n Ti 3 

m = E ( <wo = E E E = k > lo sM, 

i=l i =1 t =1 k =0 

with corresponding score function given by 

n n Ti 3 




i=l 


i =1 t =1 k =0 


I(Rit = fc) dXuk 

\tk d^ T 


The maximum likelihood estimator i/y is obtained by solving So^) = 0. 

For the missing covariate model, the observed likelihood function for 7 is 


n „ 

£ 3 ( 7 ) = 11 / Pr{Xi\Zi,Yi)dXT, 

i= 1 ■' 


with score function S 3 ( 7 ) = YJi=\ S 3 ih), where 6 ^( 7 ) = Slog / Pr{Xi\Z h Y^dXf/dj T . 
Similarly, a consistent estimator of 7 can be obtained by solving S 3 ( 7 ) = 0. 
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4.2 Estimation and Inference for the Doubly Robust Method 

Let’s denote the vector of all parameters as 0 = (/3 T , 'y T ) T ■ Our primary interest 
lies is in estimating /3. Such task can be accomplished by plugging in the estimates and 
7 in Q and solving the estimating equations for /3, that is, 

n 

Si(/3,^,7) = ySu(/3,<A,7) = 0. (10) 

1=1 

Second term of Su can be written, for X discrete, as 


E(Y?,x?\Y?,x?,Zi,Ri) {DiNi(Yi - //*)} 


^ ^ ixy {E)iNi(Yi ■ 

(wr.*D 


where the weight Wi xy is given by 


W 


ixy 


= P(YT = vTlY^Xi = Xi, Zi] (3*) x P(X r = X? 


\Y°i,X C i 


Zu 7 ). 


In the case of X continuous, the second term in Su takes the form 


E ( y~,x 7 ‘Iy-,x-,z„r ,) {DiJVify, - li,)} = f *>«, {D,Ni(Y t - M ,)} 


with conditional probability 


w 


ixy 


= P(YT = y7\Y° , X, = a*, z i5 p) X P(X™ = x;, z i5 7 ). 


The expectation in ([9]) can be cumbersome, depending on the missing data pattern. In 
such case, instead of using numerical integration, a Monte Carlo method can be applied 
to approximate the corresponding integral. 

In this work it is assumed independence working correlation and it is adopted a sand¬ 
wich standard error as given in Appendix. 


5 Simulation Study 

A small simulation study taking into account different sample sizes was conducted in 
order to quantify the bias and precision under misspecihcation of the predictive models. 
It is considered a study with = T = 3 repeated ordinal measures (with three categories) 
and two covariates (quantitative and qualitative). The true marginal model is 

logit Pr[Oit < j\X it , Z it ) = @ 0 j + f3\X it + p 2 Zu, j = 1, 2. (11) 

where Z lt is normal with unity variance and mean (0, 0.5,1) for t — 1,2, 3. 

The binary covariate X lt may be missing at some time points and is generated accord¬ 
ing to 

logit Pr{X it = l\X it , Z it ) =70 + 7 iX m _i + 72 Z it . (12) 

It is assumed <3 0 i = —0.4, /3 0 2 = 1.2, ,5] = —0.5, (3 2 = 0.5, y 0 = log{ 1), 71 = 2 and 72 = 2. 
The correlated ordinal responses were generated according to the algorithm proposed by 


8 


Touloumis (Anestis Touloumis 2013) with constant correlation between the latent vectors 


set equal to p = 0.9. 

As independent estimating equations were fitted, Ru can be defined as the indicator 
of observing both O it and X it , and it was taken 


log 


PrjRit = 1 ) 

Pr(R it = 0) 


— 0ot + ^’iI(Ri,t- 1 — 1) + ^O* t _ i + 0 sX* t _i + ipiZit, t — 2, 3, (13) 


where 0* t _ x = O i:t -i, if is observed and 0 otherwise, and X* t _ ± = A0 t _iif X it _\ 

is observed and 0 otherwise. The true values are taken as 002 = 6.6, -003 = 6, -0i = 2, 
-02 = —2, -03 = —2 and -04 = 2. ft was observed about 24% of missing observations under 
this setup. 

For comparison purposes, it was considered ordinary GEE for the complete and avail¬ 
able data, respectively, weighted GEE (WGEE), multiple imputation (MIGEE) by chained 


equations (van Buuren & Groothuis-Oudshoorn, 2011) with M = 10, and the proposed 


doubly robust version (DRGEE). In order to investigate robustness of these methods, the 
predicted models were also misspecified by omitting the covariate Ab_i from the covariate 


model (12) or the missing data model 


Results are summarized in Table [I] In each of the S = 1000 Monte Carlo replications it 
was obtained the relative bias percentage for each parameter, defined as 100 x (f3—/3) / (3, its 
standard deviation obtained through the sandwich estimator, and the coverage probability 
as a nominal 95%. 

Specially, MAR missingness impact over the response and the covariate is observed for 
all the regression parameters, the largest relative bias occur in binary covariate X. This 
comes in addition to the natural increase of parameter uncertainty. Bias in the intercept 
coefficients imply incorrect predicted probabilities for the levels of the response, whereas 
bias for parameter estimates associated with the regression covariates may erroneously at¬ 
tenuate or highlight an effect, thus leading to misinterpretations related to the importance 
of a given predictor on the longitudinal dynamics of the ordinal response. 

It can be observed that for small sample size even the GEE with for complete data 
presents a certain degree of bias. Increasing the sample size allows to clarify the perfor¬ 
mance distinctions among the compared methods. WGEE and MIGEE methods are valid 
when the model for the weight or the imputation model, respectively, are correctly spec¬ 
ified. In this case it is noted that both methods give good results for large sample sizes, 
the main distinction between them being due to the greater variability of the estimates 
for the weighted estimator. 

Doubly robust method requires the simultaneous specification of two predictive model. 
When at least one of them is correctly specified the resulting estimator is still consistent. 
Estimates are, on average, closer to those obtained with fully observed data compared to 
WGEE or MIGEE. This behavior is systematic and it can be observed to all parameters. 
By increasing the sample size the estimates from DRGEE present empirical bias in general 
smaller than their single robust competitors. This is specially true for the parameter 
associated with the incomplete binary variable X. Regarding the uncertainty of parameter 
estimates, it is noted that the variability in DRGEE is greater than multiple imputation, 
but of the same order as of the weighted method. Further, the efficiency of the doubly 
robust estimates appears relatively more sensitive to misspecihcation of the weight model 
than the covariate model. Empirical coverage rates were acceptable for correctly specified 
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Table 1: Relative bias percentage, standard deviation and empirical coverage for 1000 
simulations of incomplete covariate and response data. 


Empirical Bias Standard Deviation Empirical Coverage 



Aoi 

Ao 2 

Ai 

A2 

An 

Ao 2 

Ai 

A2 

Aoi 

A 02 

Ai 

A2 

Complete 

-4.60 

7.89 

15.07 

7.38 

0.379 

n = 50 
0.410 

0.458 

0.179 

0.94 

0.94 

0.97 

0.95 

Available 

-19.81 

15.37 

-8.36 

-11.66 

0.386 

0.420 

0.473 

0.190 

0.94 

0.92 

0.97 

0.94 

WGEE(r+) 

-12.33 

12.02 

19.96 

1.74 

0.406 

0.440 

0.513 

0.204 

0.94 

0.94 

0.97 

0.95 

WGEE(r-) 

1.50 

7.43 

1.31 

0.51 

0.413 

0.450 

0.508 

0.201 

0.95 

0.93 

0.96 

0.95 

MIGEE(z+) 

22.17 

-3.44 

-38.26 

-4.73 

0.386 

0.417 

0.480 

0.186 

0.96 

0.96 

0.97 

0.96 

MIGEE(a: _ ) 

15.81 

-0.13 

-26.18 

-2.54 

0.387 

0.417 

0.481 

0.187 

0.95 

0.96 

0.97 

0.96 

DRGEE(x+,r+) 

-8.86 

9.67 

22.20 

8.99 

0.410 

0.442 

0.513 

0.199 

0.95 

0.94 

0.97 

0.95 

DRGEE(o;-,r+) 

-7.77 

9.44 

19.33 

7.91 

0.471 

0.507 

0.579 

0.204 

0.95 

0.94 

0.97 

0.95 

DRGEE(a;+,r-) 

-7.95 

10.05 

20.58 

8.26 

0.433 

0.473 

0.540 

0.197 

0.95 

0.93 

0.96 

0.95 

DRGEE(^-,r-) 

2.86 

6.42 

-0.28 

2.99 

0.414 

0.450 
n = 150 

0.529 

0.200 

0.95 

0.93 

0.97 

0.95 

Complete 

-1.36 

1.53 

5.13 

1.99 

0.220 

0.236 

0.264 

0.103 

0.94 

0.94 

0.96 

0.95 

Available 

-17.50 

8.94 

-18.41 

-17.15 

0.224 

0.242 

0.273 

0.109 

0.91 

0.92 

0.94 

0.91 

WGEE(r+) 

-5.84 

3.41 

13.47 

1.36 

0.247 

0.267 

0.317 

0.126 

0.93 

0.93 

0.95 

0.94 

WGEE(r-) 

10.20 

-2.04 

-12.91 

-2.59 

0.253 

0.279 

0.307 

0.120 

0.94 

0.93 

0.93 

0.94 

MIGEE(z+) 

8.22 

-2.57 

-14.36 

-2.40 

0.228 

0.244 

0.283 

0.108 

0.94 

0.94 

0.96 

0.95 

MIGEE(x _ ) 

9.08 

-1.97 

-16.33 

-3.60 

0.226 

0.242 

0.280 

0.108 

0.93 

0.93 

0.95 

0.94 

DRGEE(z+,r+) 

-4.21 

2.64 

11.90 

4.12 

0.249 

0.269 

0.317 

0.117 

0.94 

0.94 

0.96 

0.95 

DRGEE(x _ ,r+) 

-5.02 

3.07 

12.26 

3.63 

0.320 

0.345 

0.396 

0.119 

0.94 

0.94 

0.95 

0.95 

DRGEE(x+, r“) 

-1.86 

1.99 

7.26 

2.95 

0.249 

0.275 

0.313 

0.113 

0.94 

0.94 

0.95 

0.94 

DRGEE(^-,r-) 

9.88 

-1.86 

-15.60 

-2.86 

0.246 

0.269 
n = 300 

0.317 

0.116 

0.94 

0.93 

0.94 

0.94 

Complete 

1.42 

0.19 

-0.04 

0.85 

0.156 

0.167 

0.187 

0.072 

0.96 

0.94 

0.94 

0.95 

Available 

-14.82 

7.67 

-23.24 

-18.14 

0.159 

0.171 

0.194 

0.077 

0.93 

0.91 

0.92 

0.87 

WGEE(r+) 

-0.24 

0.75 

5.61 

1.23 

0.182 

0.198 

0.235 

0.095 

0.95 

0.94 

0.94 

0.93 

WGEE(r-) 

15.64 

-4.50 

-22.17 

-3.89 

0.184 

0.205 

0.222 

0.088 

0.95 

0.94 

0.92 

0.94 

MIGEE(z+) 

6.23 

-1.86 

-9.95 

-1.44 

0.162 

0.174 

0.201 

0.076 

0.95 

0.94 

0.94 

0.95 

MIGEE(a: _ ) 

10.27 

-2.70 

-18.22 

-3.89 

0.160 

0.172 

0.200 

0.076 

0.94 

0.92 

0.92 

0.94 

DRGEE(i+,r+) 

-0.06 

0.83 

2.82 

1.28 

0.185 

0.202 

0.242 

0.089 

0.96 

0.95 

0.95 

0.94 

DRGEE(o;-,r+) 

-0.71 

1.08 

3.84 

1.43 

0.194 

0.213 

0.251 

0.088 

0.96 

0.95 

0.95 

0.94 

DRGEE(a:+,r-) 

2.53 

0.03 

-2.11 

0.12 

0.178 

0.198 

0.224 

0.082 

0.96 

0.96 

0.94 

0.94 

DRGEE(^-,r-) 

13.86 

-3.69 

-24.29 

-5.51 

0.176 

0.194 
n = 600 

0.228 

0.084 

0.95 

0.94 

0.93 

0.94 

Complete 

-1.01 

0.57 

1.93 

0.82 

0.110 

0.118 

0.132 

0.051 

0.95 

0.93 

0.94 

0.94 

Available 

-16.78 

7.85 

-22.20 

-18.43 

0.112 

0.121 

0.137 

0.054 

0.91 

0.88 

0.92 

0.78 

WGEE(r+) 

-1.00 

0.28 

5.43 

1.25 

0.132 

0.144 

0.170 

0.070 

0.95 

0.94 

0.95 

0.95 

WGEE(r-) 

14.87 

-4.99 

-21.53 

-3.03 

0.133 

0.149 

0.158 

0.064 

0.94 

0.94 

0.92 

0.93 

MIGEE(z+) 

1.76 

-0.66 

-3.63 

-0.39 

0.115 

0.123 

0.143 

0.054 

0.95 

0.94 

0.95 

0.94 

MIGEE(a: _ ) 

6.96 

-2.01 

-14.18 

-3.38 

0.113 

0.122 

0.141 

0.054 

0.94 

0.93 

0.93 

0.93 

DRGEE(z+,r+) 

-1.49 

0.59 

3.19 

0.98 

0.132 

0.144 

0.174 

0.062 

0.96 

0.94 

0.94 

0.94 

DRGEE(x _ ,r+) 

-1.39 

0.56 

2.96 

0.91 

0.136 

0.149 

0.178 

0.062 

0.96 

0.94 

0.96 

0.95 

DRGEE(a;+,r-) 

0.42 

0.08 

-0.12 

0.47 

0.126 

0.140 

0.159 

0.058 

0.96 

0.95 

0.95 

0.93 

DRGEE(^-,r-) 

12.24 

-3.82 

-23.18 

-5.35 

0.124 

0.138 

0.161 

0.060 

0.93 

0.93 

0.92 

0.93 


indicates correctly specified model and “ ” indicates misspecified model omitting the Xt predictor 
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Figure 1: Boxplot of the relative bias for parameter estimates under correctly specified 
models 

WGEE and MIGEE as well as for DRGEE when at least one of predictive models are 
correctly specified. 

Figure [l] shows boxplots of the percentage relative bias for the methods expected to be 
valid. The boxs represent, respectively, GEE for complete and available data, WGEE(r + ), 
MIGEE(x+), DRGEE(x + , r + ), DRGEE(x“, r + ) and DRGEE (x+,r~). As the degree of 
bias is different in the parameter estimates, for easy of visualization, they are represented 
in different scales. It can be seen that the estimates with larger variability were those 
associated with the first intercept and the covariate with incomplete values. Proposed 
method presents median relative bias close to zero and very similar to those observed 
with complete data for all sample sizes. Variability of the DR estimators is slightly larger 
than MI but of the same order as the weighted estimator. As expected, for all methods it 
can be noticed a decrease in the relative bias with increasing sample size, reflecting their 
theoretical asymptotic consistency. 

Figure [2] allows the comparison of methods incorrectly specified. The boxs repre¬ 
sent, respectively, GEE for complete and available data, WGEE(r“), MIGEE(x“), and 
DRGEE {x-,r~). 

When a key covariate is omitted from the weight model and/or the imputation one, it is 
expected all methods to be biased. This is specially true for the first intercept and for the 
incomplete binary covariate. Bias for MIGEEfT - ) seemed to get smaller than WGEE(r _ ) 
as sample sizes increases. Bias for DRGEE(a: _ , r~) was comparable to WGEE(r _ ) and 
slightly higher than multiple imputation for large sample sizes. In terms of variability of 
the estimates, the same pattern is observed as with the correctly specified models. That 
is, the multiple imputation is more efficient, followed by the doubly robust estimator and 
the weighed estimator. 
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Figure 2: Boxplot of the relative bias for the estimates with incorrectly specified models 

6 Data Analysis: Analgesia in Childbirth 

This study was conducted in Minas Gerais state, Brazil, in order to compare two 
techniques of analgesia for labor pain. There were 49 patients who were monitored during 
their entire labor period until childbirth. Pain intensity was subjectively assessed by 
the patient, and measurements of blood pressure, maternal heart rate, consumption of 
oxytocin, sedation level, signs of respiratory depression, apnea, and other variables were 
recorded. One of the techniques used was epidural analgesia (the gold standard), which 
is a local anesthetic. The other one, whose efficiency was to be compared to the gold 
standard, involved continuous intravenous infusion of remifentanil. 

The response of interest is the intensity of pain as measured by a Visual Analog Scale 
(VAS) (1: tolerable and mild pain; 2: moderate pain that causes discomfort; 3: intense and 
unbearable pain). Three measurements (0, 60,90 minutes) were selected for data analysis. 
Predictor variables considered were treatment GROUP (0: peridural; 1: remifentanil), 
AGE (in years), DU (uterine dilatation), and OXYT (consumption of oxytocin). The 
OXYT is a time-varying ordinal covariate, coded as 1, if no consumption, 2, if consumption 
equals to 10 or 30, and 3, if consumption equals or above 45. Thes eother covariates were 
chosen after a previous exploratory analysis. 

The response and oxytocin consumption was missing for 9 patients at the time 60 and 
for 18 patients at time 90. Missing is due to childbirth happened before 60 or 90 minutes. 
Therefore a MAR mechanism seems to be a reasonable assumption for this data set. The 
others covariates in the analysis were fully observed. 

For the ordinal response it was used the following proportional odds model 

logit Pr(PAIN itj < j\u it ) = /3 0j + uj t (3, j = 1,2, t = 1, 2, 3, (14) 

where u lt is the covariate vector at time t, and it is formed by TIME, GROUP, AGE, 
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Table 2: Regression Parameters for the Analgesia in Birth Data 



Available 

WGEE 


MIGEE 


DRGEE 


Parameter 

Est SE P 

Est SE 

P 

Est SE 

P 

Est SE 

P 


INTERCEPT1 

1.796 

1.308 

0.170 

1.827 

1.267 

INTERCEPT2 

3.302 

1.314 

0.012 

3.261 

1.310 

TIME 

-0.182 

0.286 

0.524 

-0.157 

0.315 

GROUP 

-1.221 

0.445 

0.006 

-1.244 

0.439 

AGE 

-0.066 

0.035 

0.057 

-0.072 

0.033 

DU 

-0.362 

0.158 

0.022 

-0.354 

0.168 

OXYT(=2) 

0.727 

0.554 

0.189 

0.712 

0.554 

OXYT(=3) 

1.285 

0.510 

0.012 

1.431 

0.503 


0.149 

1.599 

1.237 

0.196 

1.649 

1.315 

0.210 

0.013 

3.105 

1.242 

0.012 

3.069 

1.351 

0.023 

0.619 

-0.110 

0.273 

0.687 

-0.192 

0.293 

0.511 

0.005 

-1.056 

0.409 

0.010 

-1.008 

0.390 

0.010 

0.027 

-0.062 

0.031 

0.046 

-0.071 

0.035 

0.045 

0.035 

-0.372 

0.156 

0.017 

-0.340 

0.166 

0.040 

0.199 

0.954 

0.542 

0.078 

0.821 

0.531 

0.122 

0.004 

1.410 

0.488 

0.004 

1.468 

0.475 

0.002 


DU and OXIT. In response model, TIME predictor was expressed in hours rather than 
minutes. 

When using WGEE or DRGGE it is necessary to correctly model tt, in order to obtain 
consistent estimates of (3. For the missing data process Ru was defined as the indicator 
of observing both PAINu and OXIT\ t , and take the following form 

l ° S ( prSl = 0) ) = + W ‘ = 2 ' 3 ’ (15) 

where Wu includes GROUP, AGE, DU, histories of OXYT and PAIN, and the previous 
indicator of missing data. 

The distribution of the missing covariate OXYT also needs to be specified in a predic¬ 
tive model. With this aim, it was assumed a proportional odds model of the form 

logito Pr(OXYT itj < j\v it ) = 7 0i + j = 1,2, t = 2, 3 (16) 


where v lt includes main effects for GROUP, AGE, and DU. It can be observed that the 
estimate of 7 is not of interest, however it is necessary to model the missing mechanism 
related to covariate as close as possible to true in order to obtain valid estimates of (3. 
The same is true for the missing data process. All predictors in this model process were 
maintained since an overspecification is better than a underspecihcation. 

Results from four methods are shown in Tabled The first one is the usual GEE method 
using the available data; the second is the weighted method (WGEE) using model (15) 
for the weights; the third is the multiple imputation by chained equation (MIGEE) in the 
R package mice; and the fourth, labeled DRGEE, is the proposed doubly robust method 
using (15) and (16) for the weight and the covariate models, respectively. It was used an 
independent working correlation. 

TIME effect is non significant for all the four methods. All methods provide the same 
conclusion for effects of GROUP. The negative effect for GROUP means that the chance 
of women feel mild pain is lower among the group receiving the remifcntanil compared to 
the peridural group (the estimated odds is e -1 ' 008 = 0.365 in the doubly robust method). 
All methods also agree with respect to the effect of DU. That is, for each increase of 1 cm 
of uterine dilation the chance of the parturient feel mild pain decreases (in the DRGEE it 
is e ~ 0 ' 340 = 0.712, for example). It can be noticed that p -value for AGE effect goes from 
a non-significant of 0.057 in the standard GEE to a significant one in DRGEE, as well 
as for the other two missing data approaches. The conclusion is that older women have 
lower chance of experiencing mild pain than young women. For the OXIT covariate all 
methods reached the same conclusion. 
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7 Discussion 


When longitudinal ordinal data are of interest, doubly robust estimator is a nice al¬ 
ternative. Doubly robust method combines ideas from weighting and imputation and has 
been applied elsewhere for estimation of means, causal inference, and in the longitudinal 


setting for binary response data (Bang & Robins (2005), Carpenter et al. (2006), Seaman 


& Copas (2009), Chen & Zhou (2011), Li et al. (2013)). However, as far as we know, it 


has not been investigated for the longitudinal ordinal case. A doubly robust estimator is 
attractive in the sense that it needs only the correct specification of at least one of the 
models, but not necessarily both. Simulation results have indicated that, when at least the 
covariate model or missing data model is correct, the doubly robust estimators are con¬ 
sistent and present small-sample bias comparable to single robust alternatives MIGEE or 
WGEE. The proposed method presented good coverage probabilities, as well as its com¬ 
petitors but with a slight larger variance than multiple imputation. Simulation results 
also indicated that the bias of doubly robust estimators when both the covariate model 
and the missing data model are incorrect was the same magnitude of misspecihed WGEE 
or MIGEE. We hope that, in practical applications, none of the predictive models would 
be grossly misspecihed and then the proposed estimator would have a great potential of 
reducing the bias if the MAR assumption is correct. 

When the assumed independent working correlation structure differs from the true 
underlying structure, there is no price to pay in terms of the consistency and asymptotic 


normality of /3, but such a poor choice may result in loss of efficiency (Molenberghs 


& Verbeke, 2005). However, modeling of the association structure in the presence of 


missing data remains a challenge, specially with longitudinal ordinal data, because there 
is no direct way of modeling the association parameters. Future research involves the 
investigation of the impact of other association structures in the doubly robust estimates. 

In the proposed doubly robust estimator marginal means were modeled by cumu¬ 
lative logits. This implies a proportional odds model that in some cases may not be 
valid. Another future possible extension of the proposed model is, therefore, to allow 


non-proportional odds for a subset of the explanatory variables (Peterson & Harrell Jr 


1990). 
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A Appendix 

A.l Asymptotic Variance 

To state the asymptotic properties of (3, let 
Su(/3,*l>, 7 ) be the individual’s contribution to the estimating equations for /3, 

S 2 i {fi)) be the individual’s contribution to the estimating equations for if, and 
S 31 ( 7 ) be the individual’s contribution to the estimating equations for 7 . 

Define V{(3, ip, 7 ) = E {dSu({3, if, 7 )/ df3 r }, Ji 2 (/3,^,7) = £ {dS u (P, *l>, l)/d^ T }, 

/= £ {dSii(/3,t/h7)/<97 r }, / 2 W = # {<9S 2i (7/>)/<9?/> T }, J 3 (7) = ,5 {<9S 3i (7)/^7 r }, 
and QiiP,^, 7 ) = - I 12 (P 1 'ij) 1 j)If\^)S 2i (^) - I 13 (P 1 ^ 1 -Y)If 1 (-y)S 3i ^). 

Theorem 1 If either the missing data model or the covariate model is correctly specified, 
then 

n l /\p-(3 0 )^N(0, r- 1 (/3 o ,'0 o ,7o)S{r- 1 (/3 o ,'0o,7o)} T ) ) (17) 

where (3 0 is the true value of (3, and 7 0 are the probability limits of xf> and 7 , and 

s = E{Qi(Po,'l>o,'Vo)Qf(Po,'l’o,'Vo)}- 


Inferences for (3 follows by replacing the unknown quantities in (17) by its consistent 
estimators. We make use of “generalized information equality” (Pierce, 1982) that 
E{dSu(P,tl>, 7 ) /dfij T } = -E {S ll (p,'il>,~f)Sl l {'ip)}, and 


(Robins et al. 

1995), 

= -Var{S 3i ( 7 

}■ 


The matrix T is replaced by T — n 1 ELi fdSu( 6 )/dp }>, and S by X = n 1 EILi 1 QiQi 

^ ~ , - --1 




Q, = 5 M (0)-J 12 (0)J 2 ^)S 2 ^)-I 13 (0)I 3 ^‘(7)^(7), / 12 (e) = n- 1 YTi=x \ dS u {0)/dV 

/13(e) = n _1 Er=i {eS’i. t (e)/c)7 T |, I 2 ($) = n~ l EE {0S 2i ($)/0^ T }, 

h{l) = ™ _1 EE {c>5 3 i(7)/cl7 T }^ 


The proof is similar to Chen & Zhou (2011) and is omitted here. 
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